Transient Rectification of Brownian Diffusion with Asymmetric Initial Distribution 
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In an ensemble of non-interacting Brownian particles, a finite systematic average velocity may 
temporarily develop, even if it is zero initially. The effect originates from a small nonlinear correction 
to the dissipative force, causing the equation for the first moment of velocity to couple to moments 
of higher order. The effect may be relevant when a complex system dissociates in a viscous medium 
with conservation of momentum. 
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I. INTRODUCTION 

Stochastic processes with nonlinear dissipation are 
ubiquitous and challenging to describe theoretically. 
Mathematical difficulties related to the nonlinearity of 
a corresponding stochastic differential equation are only 
part of the problem. A more subtle challenge is to estab- 
lish fluctuation-dissipation relations which, in contrast 
to linear processes, cannot be phenomenologically justi- 
fied Instead, a truly dynamical approach is usually 
needed when the dissipation force and statistical prop- 
erties of the noise are deduced directly from underlying 
dynamics, rather than postulated ad hoc. Conventional 
assumptions of a phenomenological approach in the con- 
text of nonlinear response may be misleading. For in- 
stance, the assumption of Gaussian random force in the 
Langevin equation leads to the Fokker-Planck equation of 
second order, regardless of whether the dissipation force 
is linear or not. On the other hand, a kinetic approach 
leads to the second-order Fokker-Planck equation for a 
Brownian particle only in the lowest order of a pertur- 
bation technique, while in general the equation involves 
derivatives of order higher than two [E @, H, 0| • 

Nonlinear stochastic processes are usually associated 
with far-from-equilibrium dynamics. If a system is close 
to equilibrium, nonlinear dissipation usually appears as 
small corrections to the dominating linear friction and in 
many cases may be safely neglected. However, under cer- 
tain circumstances, the contribution of linear terms may 
vanish identically or be strongly reduced. Then nonlinear 
dissipative effects come into the limelight and give rise to 
a variety of new physical effects. 

An example, which has received particular attention in 
recent years, is the rectification of thermal fluctuations in 
the so-called adiabatic piston problem [5,]. The problem 
concerns Brownian motion of a piston which separates 
a gas-filled cylinder into two compartments with differ- 
ent temperatures and gas densities. If the pressure on 
both sides of the piston is the same, the linear theory 
predicts zero average velocity of the piston, while the 
correct result is that the piston acquires a systematic 
average speed in the direction of the compartment with 
higher temperature. The effect may be readily explained 
using the Langevin equation with a small nonlinear cor- 




FIG. 1: Initial velocity distribution for an ensemble of Brow- 
nian particles, discussed in the paper. The widths and heights 
of the distribution's wings are chosen so that the average ini- 
tial velocity (V) of the ensemble is zero, but the higher mo- 
ments (V n ) are finite. 



rcction, quadratic in the piston's velocity, to the dissipa- 
tive force Q . Some other effects related to the nonlinear 
dissipation are discussed in 

In the adiabatic piston problem the fluctuation- 
induced drift originates from noncquilibrium and asym- 
metry. The current point of view is that these two ingre- 
dients are necessary in general for rectification of thermal 
fluctuations, i.e. for the physical realization of Maxwell's 
demon. Asymmetry may be introduced by surroundings, 
as in the adiabatic piston problem, or by the geometry 
of the Brownian particle itself [g, In this paper, our 
concern is a transient rectification effect originating from 
asymmetric initial conditions. 



II. THE PROBLEM 

Consider an ensemble of non-interacting Brownian par- 
ticles diffusing in one dimension. The particles are identi- 
cal but may have different initial velocities. Suppose the 
distribution of initial velocities fo(V) is similar to Fig. 
1: asymmetric but in such a way that the average initial 
velocity of the ensemble is zero, 



(V(Q)) = J dVf (V)V = 0. (1) 

The question is whether (V(t)) for later time t > is 
positive, negative or zero? 
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FIG. 2: Simulation (solid) and theoretical (dashed) curves 
for the time dependence of the average velocity (V(t)) of an 
ensemble with an initial distribution similar to Fig. 1. The 
molecule-particle mass ratio parameter is A = y/m/M = 0.1. 
The widths of the distribution wings are Vi = 1/4 and V2 = 
1/2. Velocity is in units vth = y/kT/m and time is in units 
r = (A 2 7o )- 1 . 



o 
o 
<u 
> 

O) 

to 




FIG. 3: Same as for Fig. 2, but for an initial velocity dis- 
tribution with widths Vi = 1 and V2 = 2 (in units Vth) for 
the right and left wings, respectively. The corresponding en- 
semble is far from equilibrium, Vi,V2 > A. The theoretical 
(dashed) curve, given by Eq. (|16[1 overestimates the result of 
the simulation (solid curve). 



Contrary to its apparent simplicity, the question re- 
quires going beyond the standard theory of Brownian mo- 
tion based on the linear Langevin equation and the cor- 
responding second-order Fokker-Planck equation. Both 
approaches give the linear relaxation law dt(V(t)) — 
—j(V(t)), and therefore predict that if the average ve- 
locity (V(t)) is zero initially, it remains so later on. The 
prediction is incorrect as one can see from the result of 
numerical experiment presented in Fig. 2 and Fig. 3. On 
the time scale of order r = 1/7, a finite average veloc- 
ity develops in the direction corresponding to the higher, 
narrower wing of the initial distribution, the right wing 
in Fig. 1. The particles moving to the right, have lower 
average initial speeds but are more numerous and give a 
larger contribution to (V(t)) than the particles moving to 
the left. The victory of the larger team of slower runners 
does not last very long: after reaching a peak at roughly 
one half of r, the function (V(t)) decays exponentially 
with the characteristic time of order r. Yet, this tran- 
sient time may be sufficiently long to cause measurable 
physical consequences. 

The problem may be considered as an idealized model 
of the dissociation of a complex system in a vis- 
cous medium. Among possible relevant fields are the 
Coulomb fragmentation of multiply-charged clusters and 
droplets [l(| and processes involving fragmentation of 
complex molecular aggregates, such as protein-ligand dis- 
sociation [ll|. If the system is initially at rest and all 
dissociated fragments have the same mass, Eq. ([TJ is 
just the condition of conservation of total momentum. 
For a system in vacuum, the speed of the center of mass 
of fragments remains zero after dissociation. However 
if dissociation happens in a viscous medium, the aver- 
age velocity is temporarily finite, and the center of mass 
changes position even if the fragments are identical and 
have the same diffusion coefficients. 



To account for this transient rectification effect, one 
has to take into account that the equation for the first 
moment of the velocity d t (V(t)) = —j(V(t)) is closed 
only in lowest order in the small parameter A 2 = m/M, 
the mass ratio of a molecule (m) to a Brownian particle 
(M). At higher orders in A, the first moment (V(t)) is 
coupled to the moments of higher orders (V n (t)). If ini- 
tially the first moment is zero, but the higher moments 
are finite, as for the initial distribution in Fig. 1, then 
(V(t)) for t > 0. To describe the problem quantita- 
tively, one may adopt the approach based on either the 
Langevin equation for V{t) or the Fokker-Planck equa- 
tion for the distribution function f(V, t). In what follows, 
we discuss both approaches and outline details of the nu- 
merical simulations presented in Fig. 2 and Fig. 3. 

III. THEORY: LANGEVIN EQUATION 

The microscopic derivation of the Langevin equation 
beyond the lowest order in A = yj m/M was discussed 
recently in detail in [l3| . Here we outline the results and 
apply them to our problem. An appropriate perturbation 
technique is guided by anticipation that the velocity V of 
a Brownian particle is typically about A times that of a 
molecule of the surrounding bath. This suggests working 
with the scaled velocity of the particle v = X~ 1 V, which 
is expected to be of the same order as the thermal velocity 
of molecules Vth, 

v = X^V ~ vth = \r^-- (2) 
V m 

The microscopic equation of motion for the scaled veloc- 
ity v — V I A (or for the scaled momentum p = mv = 
XMV) involves the small parameter A explicitly, and 
therefore is convenient for a perturbation analysis. The 
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equation is coupled with bath degrees of freedom which 
may be "projected out" with an appropriate projection 
operator technique [H, E|- As a result, to lowest or- 
der in A, one obtains the conventional linear Langevin 
equation 



v(t) = -X 2 l0 v(t) + —F (t), 
m 



(3) 



where the zero-centered fluctuating force Fq (t) is related 
to the dissipation constant 70 through the fluctuation- 
dissipation relation 
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dt{F (0)F o {t)). 



(4) 



The linear Langevin equation ([3|) leads to the following 
equations for the velocity moments [3] 

^ = - A 2 n 70 (v n ) + A 2 n (n - 1) 7o v; 2 (v n ~ 2 ). (5) 

As discussed, these equations, obtained in lowest order in 
A , are not sufficient for our purpose: the closed equation 
for the first moment dt(v) — — X 2 jo(v) clearly cannot 
account for the behavior presented in Fig. 2. 

The next approximation for the Langevin equation in- 
volves a correction of order A and, for a homogeneous 
bath, has the form [l3| 

v(t) = -\ 2 ll v(t)-\ 4 l2 v 3 (t) + -F(t). (6) 

m 

Besides the presence of the nonlinear dissipative term 
— A 4 72 v 3 , this equation differs from the linear one (|3|) by 
a higher order correction to the linear damping, and the 
fluctuating force 



7i-7o + 0(A 2 ), F(t) = F {t) + O(A). 



(7) 



The explicit form of these corrections is not necessary for 
the purpose of this paper. The fluctuation-dissipation 
relation for the nonlinear dissipation coefficient 72 in- 
volves rather complicated correlation functions [l3j], and 
to the best of our knowledge, cannot be established phc- 
nomenologically. This is in contrast to the conventional 
fluctuation-dissipation relation (J4j) for the linear dissipa- 
tion coefficient 70 which can be obtained using the pre- 
diction of equilibrium statistics (v 2 (t)) — > kT/ra in the 
long time limit. 

Since the fluctuating force is zero-centered to any order 
in A, it follows from Eq. © that to order A 4 the first 
moment is coupled to the third one, 



-A 2 7 i (v)-X 4 l2 (v 3 ). 



(8) 



One has to substitute here (v 3 (t)) obtained in the lowest 
order in A which according to ([5]) satisfies the equation 
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Our interest is the solution of Eqs. (8j and ([9]) with the 
initial conditions 

(«(0))=0, (v 3 (0))^0. (10) 

Clearly, in this case (v(t)) ~ A 2 , so that the last term in 
the Eq. ([9]) can be neglected. Then, to order A 2 , the third 
moment decays exponentially (v 3 (t)) = (i; 3 (0))e~ 3A 7ot . 
Substituting this into Eq. ([5]) and recalling that 71 = 
70 + 0(A 2 ), one obtains 

(«(<)) = -X 2 -^- (v 3 (0)) e - x2 ^ ot (l - e- 2A ^ *). (11) 

270 

Recall also that v is the scaled velocity, v — V/X. For 
the true velocity V the result formally does not involve 
the small factor A 2 , 

(V(t)) = -— (V 3 {0)} e-^^il - e- 2A2 T 0t ). (12) 
27o 

However, one should keep in mind that the whole proce- 
dure applied above implies that V ~ A v t h ■ This puts a 
constraint on the the width A of the initial distribution 

/ (n 



A < Xv 



in ■ 



(13) 



-3A 2 7 (« 3 )+6A 2 7oVh»- 



(9) 



Under this constraint (V 3 (0)) is small and cannot exceed 
order A 3 v 3 h . 

For a far-from-equilibrium ensemble the above theory, 
strictly speaking, is not applicable. Yet, as one observes 
from Fig. 3, Eq. (TT2"|) predicts qualitatively correct be- 
havior also for a "hot" initial distribution with A ~ vth- 
In these cases the first moment given by Eq. (fT2")) is not 
small, (V(t)) - A . 

According to the result (fT2")l . the first moment (V(t)) 
reaches the maximum at time to = (ln3/2)r « 0.55 r 
where r = X^ 2 ^ 1 , which is seen in Fig. 2 to be in agree- 
ment with numerical simulation. To make more qualita- 
tive predictions, one needs an explicit expression for the 
ratio of the dissipative coefficients 72/701 which is the 
prefactor in Eq. (TT2")) . Since a general result for this ratio 
is unknown, in the rest of the paper we discuss a spe- 
cific model of Brownian motion - the Rayleigh model - 
for which our numerical experiment is carried out, and 
for which analytical results are available. 

In the original Rayleigh model [1, H, [1] ) a heavy Brown- 
ian particle moves in one dimension interacting with bath 
molecules through instantaneous elastic collisions, while 
molecules do not interact with one another at all. For 
this model the Fokker-Planck equation for the distribu- 
tion function f(V,t) can be readily obtained, as will be 
discussed in the next section. However, due to the singu- 
lar character of the hard- wall potential, the derivation of 
a nonlinear Langevin equation for the original Rayleigh 
model is not quite straightforward. One may instead to 
work with a generalized Rayleigh model where the par- 
ticle interacts with molecules through a continuous re- 
pulsive potential. For a low density of bath molecules 
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(when multiple collision are negligible) and for the time 
scale longer than the collision time r c , the original and 
generalized models are expected to give the same results. 
Using the generalized Rayleigh model, one obtains the 
following explicit expressions for the dissipative coeffi- 
cients 1131 



7o 



nSvth, 72 



3V2vr 



:TlSv. 



III 



(14) 



Here n is the concentration of molecules, S is the par- 
ticle's cross-section, and Vth = \J kT /m is the thermal 
velocity of molecules in the bath. It is tempting to as- 
sume that the relation 



in 



72 = 1 _2 = 

7o 6 Vth 6fcT' 



(15) 



which follows from (fT4"|) , is in fact general but we leave 
this conjecture for further studies. Substituting (|15|) into 
Eq. (p~2|) . one finally obtains 



(t/ 3 (0))e- A2 ^ t (l-e- 2A270 *). (16) 



(V(t)) 



12 kT 



Subsequently, the average displacement of the ensemble 
is 

(X) = / dt (V(t)) = - <U 3 (0)>. (17) 

Jo 18 7oA 2 «4 

The result OH for (V(t)), presented in Fig. 2 by 
dashed lines, is in good agreement with numerical simu- 
lation as long as the constraint (|13p on the initial distri- 
bution is satisfied. Before discussing details of the simu- 
lation, let us derive the results using the language of the 
Fokker-Planck equation. 



IV. THEORY: FOKKER-PLANCK EQUATION 

For the original Rayleigh model, which involves only 
binary-particle molecule collisions, the Fokker-Planck 
equation can be readily obtained using the Kramers- 
Moyal expansion of the master equation 0, H, 0, Q . To 
order A 2 , the equation has a familiar form 



df(v,t) _ x2 



at 



= X^ D 2 f(v,t), 



(18) 



where the second order differential operator D 2 reads 



8 



: v + v th' 



d 2 



(19) 



and 70 is given by (fbl"|) . This equation corresponds to 
the linear Langevin equation ((3J) and produces Eq. © 
for the moments (v n (t)) to order A 2 . The equation of 
order A 4 has the form 0, H[ 

- A 2 7o D 2 f(v, t) + A 4 7 o D A f(v, t). (20) 



where the forth-order differential operator D4 reads 



Da 



dv 6 
3 jP_ 

"2 dv 2 '' 



u ih 



dv 3 



'dv 2 

4 4 a 4 



For the first moment, Eq. (|20[) gives the following equa- 
tion 

±{v) = -A 2 70 (l - A 2 ) (v) - i A 4 7 o^ 2 (« 3 >- (22) 

Recalling the relations ([7]) and (115|) . one observes that 
the above equation is equivalent to Eq. 1(5)) derived from 
the nonlinear Langevin equation. Therefore, the Fokker- 
Planck equation ([20|) gives the same results as the non- 
linear Langevin equation ([6|). Note, however, that the 
Langevin equation © is derived directly from the Liou- 
ville equation [l3[ and is more general than the Fokker- 
Planck equation (f2T)]) . which is obtained under the as- 
sumption of binary particle- molecule collisions. 



V. SIMULATION 

In our molecular dynamics simulation, we use the gen- 
eralized Rayleigh model in which the Brownian parti- 
cle moves in one dimension interacting with molecules 
through a finite-range repulsive parabolic potential, while 
molecules do not interact with one another. In this 
model, discussed in detail in the particle- molecule 
collision time t c is finite and does not depend on the ve- 
locity of the molecule. A characteristic parameter of the 
model is N — nSv t hT c , which is an average number of 
molecules simultaneously interacting with the particle. 
In simulation, the linear molecular density nS was cho- 
sen to make N of order 1. In this case, multiple particle- 
molecule collisions are rare, and one can expect that the 
result should be close to that for the original Rayleigh 
model with instantaneous binary collisions. 

To mimic unbounded diffusion of a particle, we have 
used two sources of molecules located far from the par- 
ticle that generate a bath with a Maxwellian velocity 
distribution and a constant density. The first condition 
is easily accommodated by selecting incoming molecule 
velocities from the Boltzmann distribution 



iSv 



v th \/2ir 



exp 



(23) 



while controlling the rate of molecule generation with a 
Poisson process is one possibility that is consistent with 
the second condition. With such a velocity distribu- 
tion, the total flux at each source is <!> = L°° <j>(v)dv = 
nSvth/VZft- The Poisson distribution for the period be- 
tween molecule injections is then P(ri n ) = exp(— $i), 
which will maintain an average linear density of nS 
around the particle. 
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An ensemble of particles is emulated by performing 
multiple runs, resetting the system between each run 
with the new particle initial conditions selected from the 
appropriate distribution functions, and averaging the re- 
sults of all runs together. For a symmetric velocity dis- 
tribution function fo(v), the simulation reproduced fa- 
miliar results of linear Brownian motion including the 
exponential decay of the velocity correlation function on 
a time scale t > t c and deviation from exponential form 
for t < r ff , which is in agreement with the theory devel- 
oped in [la ]. 

Consider now an asymmetric initial distribution such 
as that shown in Fig. 1. Let x — V/vth be the dimension- 
less velocity of the particle. Also let x\ , X2 be the widths 
and Ci , C2 be the heights of the right and left wings of 
the distribution fo(x), respectively. The conditions of 
normalization J dxfo(x) — 1 and of zero first moment 
J dxfo(x)x — give 



C\X\ + C 2 X 2 — 1, C\X l 



C 2 x\ 







and therefore, 



ci = 



X 2 1 



C2 



Xl 1 



(24) 



(25) 



Xl Xi + X 2 X 2 Xl + x 2 

The theoretical prediction is given by Eq. (fT2"| . 

(x(i))=-l(x 3 (0))e- t /-(l- e - 2t /-), (26) 

where r = (A 2 7o)~ 1 , and the initial third moment, ac- 
cording to (|25"|) . equals 



(x 3 (0)) = ^(xi-x 2 ). 



(27) 



to satisfy in simulation. Since (x(t)) ~ (a; 3 (0)) < A 3 , 
one needs a very large number of runs (larger than A~ 6 ) 
to average out fluctuations and find the function (x(t)) 
with reasonable precision. On the other hand, a strongly 
non-equilibrium ensemble with the initial distribution 
widths xi,X2 ~ 1 is easier to simulate since in this case 
(x(t)} ~ 1, which requires a relatively small number of 
runs. 

The simulation has been performed for A = 0.1, N = 
nSvthTc = 1, time step At = 0.1 t c , and various parame- 
ters of the initial two- wing distribution fo(x) in Fig. 1. 
Time in Figs. 2 and 3 is given in units of velocity cor- 
relation time r = (A 2 7o)~ 1 which, according to (Ti~4"|) . is 
related to the collision time r c by t c /t = (8/v / 27r)A 2 iV. 

Fig. 2 corresponds to the initial velocity distribution 
fo(x) with left and right maximum velocities xi = 1/4 
and x 2 = 1/2, respectively. This is a close-to-equilibrium 
ensemble, xi,x 2 ~ A. For this case, Eqs. and (|2"T|) 
give ci = 8/3, c 2 = 2/3, and (ir 3 (0)) = -1/128. As dis- 
cussed above, this case requires a large number of runs to 
minimize relative fluctuations. The presented plot (solid 
line) is the average over about 5 x 10 7 runs. Despite still 
visible fluctuations, the data and theoretical prediction 
(jTHJ) are clearly in good agreement. 

Fig. 3 corresponds to the distribution with maximum 
velocities xi — 1 and x 2 — 2. In this case, ci = 2/3, 
c 2 = 1/6, and (x 3 (0)) = —0.5. The corresponding en- 
semble includes "hot" Brownian particles with initial ve- 
locities x > A (V > Xv t h), so that the major assumption 
of the theory is not satisfied. It is not surprising then 
that in this case the theoretical prediction (|26p distinctly 
overestimates the simulation curve. Qualitative theory 
for a strongly non-equilibrium ensemble remains a chal- 
lenge. 



Recall that the theory outlined in previous sections 
applies under the close-to equilibrium constraint (|13[) . 
which requires that xi and x 2 must be of order A or 
less. Note that for small A, this condition is not easy 
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